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ABSTRACT 


In  the  numerical  modeling  of  the  crack  propagation  in  dynamic  fracture  using  stationary  elements, 
a  discrete  and  sudden  release  of  node  at  the  crack  tip  creates  spurious  oscillation  in  the  kinetic  and 
strain  energy  values.  In  order  to  reduce  the  oscillation,  a  moving  node  element  was  utilized.  This 
element  can  model  a  continuous  crack  tip  movement  more  closely.  The  moving  node  element  is 
compatible  with  surrounding  regular  isoparametric  elements  and  no  remeshing  is  required  during  the 
crack  propagation.  In  addition,  two  different  central  difference  schemes  were  compared,  and  their 
results  were  almost  the  same. 
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I .  INTRODUCTION 


Inherent  flaws  in  engineering  materials  can  lead  to  a 
catastrophic  failure  due  to  unstable  crack  propagation.  By 
eliminating  the  conditions  and/or  manufacturing  defects 
(i.e.,  overloading,  fabrication)  which  can  lead  to  crack 
initiation,  catastrophic  failure  can  be  prevented.  In  many 
structural  components  absolute  prevention  can  not  always  be 
guaranteed.  For  such  structures,  catastrophic  failure  can  be 
minimized  by  a  crack  arrest  system.  To  ensure  the  proper 
incorporation  of  a  crack  arrest  system,  the  effects  of  a 
propagating  crack  must  be  understood.  These  effects  can  be 
understood  through  the  study  of  dynamic  fracture  mechanics. 

Dynamic  fracture  mechanics  can  be  applied  to  any  problem 
involving  a  body  containing  a  crack  in  which  inertia  forces 
play  an  important  role.  In  practice  two  kinds  of  dynamic 
fracture  problems  have  received  the  most  attention  [Ref.  1] . 
These  are: 

1.  bodies  with  stationary  cracks  that  are  subjected  to  a 
rapidly  varying  load,  and 

2.  bodies  under  fixed  or  slow  varying  loading  that 
contain  rapidly  moving  cracks. 
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Dynamic  crack  propagation  can  be  divided  into  two 
categories:  crack  initiation  and  crack  propagation.  A  third 
category  sometimes  used  is  crack  arrest.  There  are  different 
opinions  concerning  crack  arrest.  Kanninen  [Ref.  2]  asserted 
that  crack  arrest  is  not  a  separate  category  of  dynamic  crack 
propagation,  but  is  the  termination  phase  of  crack 
propagation . 

Dynamic  fracture  mechanics  problems  have  focused  mainly  on 
bodies  that  contain  a  rapidly  moving  crack.  As  such,  several 
numerical  dynamic  fracture  analyses  for  the  opening  mode  crack 
propagation  have  been  studied.  Kobayashi,  et  al .  [Ref.  3] 
analyzed  two  fracturing  Homalite-100  plates  by  dynamic  finite 
element  and  dynamic  photoelastic  analyses.  These  analyses  used 
the  process  of  discrete  crack- tip  advances.  The  restraining 
nodal  force  was  suddenly  released  when  the  crack- tip  reached 
the  next  adjacent  node.  As  indicated  by  Malluck  and  King  [Ref. 
4],  the  above  procedure  has  inherent  problems:  the  sudden 
release  of  a  node  can  induce  unwanted  high-frequency  of 
motions;  and  the  crack  tip  location  within  the  nodal  spacing 
can  not  be  determined.  To  reduce  this  problem  Malluck  and  King 
incorporated  a  mechanism  for  energy  release  in  their  finite 
element  analysis.  The  nodal  reaction  of  the  crack  tip  was 
gradually  reduced  as  the  crack  tip  propagated  (in  the 
continuum  view)  from  one  node  to  the  next  node. 

Kobayashi,  et  al.  (Ref.  5]  were  aware  of  the  considerable 
oscillations  in  the  calculated  dynamic  energy  release  rate. 
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They  believed  the  oscillation  was  due  to  the  sudden  release  of 
the  crack  tip  (finite  element  nodes) .  To  reduce  these  effects 
a  nodal  force  release  mechanism  was  incorporated  to  depict  a 
more  gradual  transit  of  the  crack  tip  between  adjacent  nodes. 
However,  the  assumed  nodal  release  force  is  a  somehow 
arbitrary  choice. 

As  indicated  in  [Ref.  6],  computational  procedures  for 
crack  propagation  problems  can  be  grouped  into  two  types: 
stationary  mesh  procedure  and  moving  mesh  procedure.  The  above 
analyses  may  be  categorized  as  a  stationary  mesh  procedure.  To 
predict  the  continuing  propagation  of  a  crack  in  a  discrete 
model  closely,  moving  mesh  procedures  were  introduced. 
Nishioka  and  Atluri  [Ref.  7]  analyzed  the  propagating  crac“L 
using  a  moving  singular  element  in  which  the  geometry  of  the 
element  was  altered  as  the  crack  propagated;  thereby, 
requiring  remeshing  of  the  elements  after  an  elapse  period  of 
time.  If  the  element  is  distorted  severely  the  accuracy  may  be 
degraded  as  indicated  in  [Ref.  8];  therefore,  it  was  critical 
that  remeshing  occurred  at  the  appropriate  time.  Kwon  and  Akin 
[Ref.  9]  developed  a  procedure  for  modeling  the  crack 
propagation  problem  using  a  node  moving  along  the  edge  of  the 
element.  In  this  procedure  the  element  is  not  distorted; 
therefore,  remeshing  is  not  required.  A  modification  of  the 
latter  procedure  was  incorporated  into  this  study. 

The  main  objective  of  this  study  is  to  examine  different 
numerical  modeling  techniques  for  studying  dynamic  cracJi 
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propagation.  The  main  emphasis  is  placed  on  reducing  the 
spurious  oscillation  in  the  calculated  energy  terms  which 
resulted  from  the  stationary  node  procedure.  Figure  la.  An 
eight  noded  regular  element  (Figure  lb)  was  used  in  the 
procedure.  To  reduce  the  spurious  oscillation  a  moving  node 
procedure  (Figure  2a)  using  a  moving  node  (Figure  2b)  was 
incorporated  into  the  finite  element  model. 

The  study  focused  on  the  analysis  of  a  rapidly  propagating 
crack  of  opening  mode  in  a  linearly  elastic  isotropic  body. 
Inertia  was  formulated  into  a  diagonal  mass  matrix  and  the 
numerical  time  integration  was  accomplished  using  the  two 
different  central  difference  method.  Crack  propagation  was 
simulated  by  the  sequential  release  (at  prescribe  time 
intervals)  of  the  nodes  along  the  edge  of  the  finite  element 
model.  Both  stationary  and  moving  node  techniques  were  used  to 
model  the  crack  tip  movement.  Crack  tip  stress  singularity  was 
not  represented  in  the  model.  From  these  method,  crack  opening 
displacement,  work,  strain  energy,  and  kinetic  energy  were 
calculated  and  a  comparative  analysis  was  conducted  from  the 
results . 
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Figure  1.  (a)  Modeling  of  a  propagating  crack  using  a 

stationary  node  element,  (b)  An  eight  noded  regular 
element . 


Figure  2.  (a)  Modeling  of  a  crack  propagating  using  a 

moving  node  element,  (b)  An  eight  noded  moving  node 
element . 


II.  MATHEMATICAL  DERIVATION 


A.  GOVERNING  EQUATION 

In  an  elastodynamic  problem  of  a  continuous  isotropic 
body,  certain  governing  equations  must  be  satisfied  at  all 
interior  points  of  the  body.  For  two-dimensional  problems,  the 
equations  of  motion  are: 


X  ^ 

dx  dy  ^  dt 


^F-p 


d^u 

d^c 


=  0 


dv’  „  ^  r\ 

_ -D - =0 

dx  dy  ^  dt  ^  dt^ 


(1) 


where  <jy,  and  r^y  are  the  normal  stresses  in  the  x  and  y 
direction  and  the  shear  stress,  respectively.  and  Fy  are 
the  body  forces  in  the  x  and  y  direction,  respectively,  p  is 
the  density  of  the  material,  fi  is  the  damping  coefficient  for 
the  viscous  damping,  u  and  v  are  the  displacements  in  the  x 
and  y  direction,  respectively.  And  t  denotes  the  time. 

The  constitutive  equation  (stress-strain  relations)  is 

{a}  =  [D]{€}  (2) 

where  {a}  represents  the  stress  vector  and  {e}  the  strain 
vector,  as  shown  below: 
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(3) 


{O}  =  [O,  Oy  V)’’ 

{€}  =  [€^  e, 


and  [D]  is  a  symmetric  3x3  material  property  matrix.  For  an 
isotropic  material,  [D]  is 


[D] 


1  V 


E 

(1-V-) 


V  1 
0  0 


0 

0 


1-v 


2 


(4) 


for  the  plane  stress  condition  and 


[Z?]=. 


E(l-v) 


(l*v) {l-2v) 


1 

V 


1-v 

0 


1-v 

1 


0 


0 

0 

l-2v 


2  (1-v) 


(5) 


for  the  plane  strain  condition. 

The  kinematic  equation  (strain- displacement  relation)  is 


•y 

tYxyJ 


du 
dx 
dv 
dy 
du  ^  dv\ 
dy 


(6) 
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B.  BOUNDARY  CONDITIONS 


In  addition  to  satisfying  the  90verning  equations,  the 
prescribed  conditions  on  stress  and/or  displacement  must  be 
satisfied  on  the  boundary  surface  of  the  body.  These 
conditions  are  classified  as  below.  The  essential  boundary 
conditions  are 


u=u 


v=v 


(7) 


and  the  natural  boundary  conditions  are 

xy^x"*^  ^  y^y 

where  <i>y.  and  <py  are  the  traction  in  the  x  and  y  direction, 
respectively.  And  n^^  and  n^  are  the  components  of  the  outward 
directed  unit  normal  vector  on  the  boundary  surface  in  the  x 
and  y  direction,  respectively.  Boundary  value  problems,  in 
general,  consist  of  a  combination  of  the  two  types  described 
above . 

To  transform  the  differential  equations  to  matrix 
equations,  apply  the  Galerkin  method  of  weighted  residual  to 
the  equations  of  motion  (1) : 
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(9) 


dx  dy  dc 


jo),  dQ  =0 


dy  ^dy  ■'‘ 


(j^dQ  =0 


where  and  u>2  are  test  (weighing)  functions  and  is  the 
domain  of  the  given  problem.  Because  Galerkin's  method  yields 
a  non- symmetric  matrix  equation,  apply  the  weak  formulation  to 
equation  (9)  to  get  a  symmetric  matrix.  The  results  yield  a 
symmetric  equation  because  the  governing  equation  is  self- 
adjoint.  Integrating  by  parts  the  stress  terms  of  equation  (9) 
yields : 


where  F  is  the  boundary  surface  of  the  domain  n.  Rewriting 
equation  (lO)  into  matrix  formr 


(11) 
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where  the  first  term  becomes  the  stiffness  matrix,  the  second 
term  becomes  mass  matrix,  the  third  term  becomes  the  damping 
matrix,  the  fourth  term  becomes  the  load  vector,  and  the  fifth 
term  becomes  the  body  force  vector.  Substituti.ng  the 
constitutive  equation  (2)  into  (11)  gives 


0 

<i)o 


dx 

0 

dy 

[D]< 

- 

0), 

0  ' 

d^u 

dt^ 

(0, 

0  ' 

(au) 

dc 

fa' 

0 

d(i). 

d<i}^ 

s 

-+p 

0 

<02 

d^v 

0 

<*>2 

dv' 

dy 

dx 

Yxy, 

dt^ 

[dc 

>)  do)  = 


(12) 


Substituting  the  kinematic  equation  (6)  into  (12)  gives 


/. 


f  6<*), 

dx 


0 

d<t>2 

~dy 


3<i>. 

“aF 

a<j2 

dx 


CO] 


du 
dx 
dv 
dy 
du  ^  dv\ 
dy  dx\ 


0), 


a^Li' 

du 

f«,  0  1 

at* 

1 

l"x  o] 

dt 

3 

o 

d^v 

at^J 

3 

o 

dv 

dc 

0 


)  dQ  = 


(13) 


The  matrices  and  the  load  vector  are  computed  on  the  element 
level  by  discretization  of  the  domain  £2  and  the  boundary 
surface  F; 


r=2:r* 


(14) 


11 


where  Qg  and  Tg  are  the  element  domain  and  element  boundary 
surface,  respectively. 


C.  MASS  MATRIX 

The  mass  matrix  comes  from 


f 

“1 

0  ' 

Jq« 

0 

“2 

d^u 


dt^ 

d^v 


at" 


dQ 


(15) 


In  this  study  eight  noded  isoparametric  elements  were  used  for 
the  finite  element  mesh.  The  finite  element  derivation  for 
the  eight  noded  element  is  shown  in  detail  below.  The  element 
is  shown  in  Figure  lb. 

From  variable  interpolation, 

U=Hj_U^^H2  U2  U3  +H5  U5  +//s  Ug  U7  +//g  Ug 

( 16 ) 

v=H^  V3  +H2  V2  V3  Vg  +//5  V5  Vg  V,  +i/8  Vg 


where  u^  and  v^  are  the  displacements  in  the  x  and  y  direction 
of  the  respective  node  and  are  the  shape  functions  at  the 
respective  node. 

For  the  regular  eight  noded  isoparametric  element 
[Ref.  10]  they  are 
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//,=—  (1-s)  (1-c) 

’  4 

H,  =  —  (1-s)  (1-C)  (s-c-1) 

^  4 

//,  =  —  (1-s)  (1-c)  (1-c)  (s+c-1) 

’  4 

H=—  (l-s)  (1-C)  (  c-s-1) 

(17) 

^5  =  ^  (1-C)  (1-S‘) 

W,=  —  (l-s)  (1-c^) 

®  2 

//,  =  —  (1-C)  (1-s^) 

2 

//.  =  —  (1-S)  (1-cO 
*  2 


and  for  the  eight  noded  isoparametric  moving  node  element 
[Ref.  9]  they  are 


=  ^  (1-s)  (1-t) 

//2  =  4  (l^S)  (l-t) 

=  ^  (1-s)  (l>t) 

//,  = - -  (1-S‘)  (1-t) 

2 (1-So) 

//6  =  -|  (1+s)  (l-t^) 

/£,  =  -|  (1  +  t)  (l-S^) 

/fa  =  -|  (1-s)  (1-t^) 


(18) 


Since  the  shape  functions  are  independent  the  time,  the 
acceleration  terms  are 
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u. 

<y  u 

dt^  0  W2  0  w,  0  0  H5  0  0  H,  0  //,  0  I  : 

d^v  0  0  Hx  0  H,  0  0  //,  0  0  //,  0  H,  I  : 

dC^  u. 


(19) 


Since  the  shape  functions  are  expressed  in  the  natural 
coordinates,  it  is  necessary  to  carry  out  the  integration  in 
equation  (21)  in  the  natural  coordinate,  using  the 
relationship 

dQ=dxdy=det  [J]  dsdt  (22) 
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where  [J]  is  Che  Jacobian  macrix  given  as 

dx  By 
Bs  Bs 
Bx  By 
Be  Be 


SubstiCuting  equation  (22)  into  (21)  yields 

u, 

V; 

U5 
V; 

dec  [  J]  dsde  ' 

Since  the  shape  functions  are  an  explicit  function  of  s 
and  t,  a  numerical  method  is  used  to  evaluate  equation  (24) . 
The  two-point  Gaussian  quadrature  formula  is  used  to  carry  out 
Che  integration.  The  result  is  a  symmetric  consistent  mass 
matrix. 

For  this  study  the  consistent  mass  matrix  was  diagonalized 
by  following  procedure  [Ref.  11] : 

a.  Calculate  the  diagonal  coefficients  of  the  consistent 
mass  maCrix. 


0 

0 

.. 

0 

0 

0 

.. 

0 

0 

0 

.  .  H2H. 

0 

0 

0 

.. 

0 

0 

0 

.. 

0 

0 

0 

.. 

0 

0 

0 

.. 

0 

0 

0 

0 

^9^8. 

(23) 
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b.  Determine  the  total  mass  of  the  element,  m. 

c.  Add  the  diagonal  coefficients,  associated  with 

translation  only  to  obtain  a  number,  e. 

d.  Multiply  the  diagonal  coefficients,  by  m/e. 

e.  Set  all  the  off-diagonal  terms  to  zeros. 


D.  DAMPING  MATRIX 

The  damping  matrix  is  given  by 


0 


du 


dt 

dv 


dt 


dQ 


(25) 


Since  the  shape  functions  are  independent  of  time. 


du 

dt 

0 

H2  0 

0 

0  //g  0  Hg  0  H,  0 

dv 

[at] 

0 

0  H2 

0 

^3 

0 

0 

0 

0 

//g  0 

0  H, 


(26) 


Substituting  equations  (20),  (22)  and  (26)  into  (25)  and 

applying  the  two-point  Gaussian  formula  gives  a  symmetric 
consistent  damping  matrix. 

E.  STIFFNESS  MATRIX 

The  stiffness  matrix  is  given  as 
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d(j>. 

dx 


0 


8(i), 

dy 


8(0, 


dy 

8(0- 


8x 


[D] 


< 


du 

dx 

dv 

dy 


}dQ 


du  ^  dv\ 

dy  8xJ 


(27) 


Rewriting  the  kinematic  equation  in  terms  of  nodal 


displacements  gives 


du 

dx 


I  ^  \ 

dy 

du  ^  ^ 
dy  dx 


0 


dy 


ib 

8x 


'  dH^ 

0 

8h, 

0  .  . 

r\ 

dH^ 

r\ 

dHg 

0 

u._ 

dx 

dx 

.  .  u 

dx 

U 

dx 

V. 

0 

dH^ 

dy 

0 

dH^ 

dy 

.  .  0 

dH^ 

~dy 

0 

dH, 

dy 

‘ 

dH^ 

dHj^ 

dH. 

dH^ 

dhj 

dH, 

8H3 

dh’g 

~dy 

8x 

dx 

dx 

dy 

dx 

(28) 


where 


dx  “  ds 


A 

dt 

dH, 


dt 


(29) 


and  A^j  are  the  components  of  the  inverse  of  the  Jacobian 
matrix.  The  matrix  in  e(3uation  (28)  is  called  the  [B]  matrix 
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and  the  vector  is  called  the  displacement  vector  (d}  .  From 
Galerkin's  method  on  the  first  matrix  in  equation  (27)  gives 


3(i), 

dx 

0 


0 

dui, 

■37 


3(0. 

dy 

3(02 

dx 


IBJ  ■ 


(30) 


Substituting  equations  (29)  and  (30)  into  (28)  gives 


f  [BldQid] 

J  Q- 


(31) 


The  stiffness  matrix  is  given  as 

[B]^lD][B]dQ  (32) 

J  Q« 

Substituting  equation  (22)  into  (32)  gives 

[iC^]  [fl]  ^[£?]  [Bidet  [^]  dsdc  (33) 

The  two -point  Gaussian  (Quadrature  formula  is  used  to  carry  out 
the  integration, 

P.  BODY  FORCE  VECTOR 

The  body  force  vector  is  given  as 
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Substituting  equations  (20)  and  (22)  into  (34)  gives 


iH,  0  0 

-iJ-i  0  0  H, 


0  He 

H,  0 


0 

^9 


det  [  J]  dsdc 


(35) 


The  two- point  Gaussian  quadrature  formula  is  used  to  carry  out 
the  integration. 


G .  LOAD  VECTOR 

The  load  vector  is  given  by 


(36) 


where  T  is  the  boundary  surface  of  the  domain  0.  If  there  is 
no  traction 


(37) 


If  there  is  traction,  substituting  equation  (20)  into  (36) 
gives 


(P)=/-  pi  °  ^  ° 

J4-7-3  Q  TT^  0  K  0  Ju\  |4>J  2 


—  (h)  ds 


where  {P}  is  Che  external  load  vector  along  edge  4-7-3  (Figure 
3)  ,  h  is  the  element  thickness  and  overbars  indicate  the  shape 
functions  are  evaluated  at  c=l  (the  natural  coordinate) . 
Integrating  equation  (38)  yields 


'"thi  i  II 


Fraction  of  the  total  force  is  shown  in  Figure  3. 


Ty-  T;^  -  2/3 

0..-  0„-  1/6 


Figure  3 .  Fraction  of  total 
force  on  the  element 
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From  the  procedural  derivation  above,  the  equation  of 
motion  is  finally  rewritten  as 

[M]  {d}-[/n  {d}  =  {;?}  (40) 

where  {r}={p}+{f}.  This  equation  applies  to  any  dynamic 
problem  in  which  the  finite  element  method  is  used. 

H.  TRANSIENT  ANALYSIS 

In  dynamic  problems  the  displacements,  velocities, 
accelerations,  strains,  stresses,  and  loads  are  time 
dependent.  Therefore,  a  time  integration  scheme  is  required. 
The  equation  of  motion  at  time  t  is 

[M]  {d]^={R}^  (41) 

This  study  examined  two  different  forms  of  the  central 
difference  method: 

The  first  method  is  as  summarized  below  [Ref.  12]; 

A.  Initial  Calculations 

I.  Solve  for  [M] ,  [C] ,  and  [K] 

2.  {d}''  =  [M]-i[{R}0  -  [C]{d}0  -  [K]{d}°]. 

where  {d}°  and  {d}°  are  the  initial  conditions. 

3.  {d}-^*^  =  {d}0  -  (At){d}0  +  a3{d}0. 

4.  [M]  =  ao[M]  +  ai[C] . 
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B.  For  each  time  step 

1.  {R}  =  {R}’'  -  (  [K]  -  a2[M]){u}''  -  (ao[M]  -  a^  [C]  ) 

2.  =  [M] 

3.  {d}*^  =  ao[{d}t"^‘^  -  2{d}'^  .  {d}"-^*^]  . 

4.  {d}*^  =  -  {d}'^-^'^] 

where  ag  =  1/At^;  a^  =  l/2At;  a2  =  2ao;  a3  =  l/a^ 

The  second  form  "summed  form"  [Ref.  13]  is  as  follows: 
A.  For  each  time  step 

1.  {d}^'^"'^  =  [MI'Mr^'^*^  -  [K]  {d}^'^"’M  . 

2.  ^  |^|At  +  l/2  ^  + 

3.  {d}^‘'-"2  ,  {d}^t  +  l  ^  ^tAt-H3/2|^}At  +  3/2^ 


The  explicit  central  difference  scheme  is  conditionally 
stable.  The  time  step  At  must  be  controlled  from  the  smallest 
period  of  the  system.  As  indicated  in  [Ref.  14],  a  safe 
limiting  value  is  given  by 


At=0.45L" 


p  (l-»-v)  (l-2v) 
E{l-v) 


1/2 


(42) 


where  L"  is  the  smallest  distance  between  adjacent  nodes  of 
any  element  used  .'>rid  E  and  «'  are  the  elastic  modulus  and 
poisson's  ratio,  respectively.  This  safe  limiting  value  was 
used  as  a  criteria  to  ensure  stability. 
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I.  DYNAMIC  FRACTURE  ANALYSIS 

From  elasticity  the  stresses  near  a  crack  tip  for  the 
opening  mode  are: 


-^cos— (l  -sin  —  sin  — j 

2WI  2  2  I 


v'2n 

K 
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K  6/,  .  0  30 

—  cos—  1  +sin—  sin  — 

Tir  2\  2  2 


K  I  ■  6  6  30 

— - — sin  — cos  — cos  — 

v/ItT?  ',222 


(43) 


where  and  Cy  are  the  tensile  stresses  in  the  x  and  y 
direction,  respectively,  is  the  shear  stress,  r  is  the 
radial  distance  from  the  crack  tip,  0  is  the  angle  from  the 
crack  plane,  and  K  is  the  stress  intensity  factor.  The  stress 
intensity  factor  provides  a  convenient  parameter  of  the  stress 
distribution  around  a  flaw.  If  K  is  equal  to  or  exceeds  a 
critical  value  (fracture  toughness,  K^.)  the  crack  becomes 
unstable.  is  a  material  property.  [Ref.  15] 

In  the  dynamic  generalization  of  linear  elastic  fracture 
mechanics  (LEFM)  ,  has  two  parts.  For  the  initiation  of 
crack  growth 

K{a,o,t)  =  K^{o)  (44) 


where  a  represents  the  crack  length,  a  is  the  applied  stress, 
t  is  time,  and  a  is  the  loading  rate.  Similarly,  for  a 
propagation  crack 
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Kia,  a .  t)  =  K^iA) 


(45) 


where  a  denotes  the  crack  speed.  [Ref.  2] 

An  alternative  criterion  for  determining  the  motion  of  a 
crack  is, 

G(a,o,  t)  =  i?(i)  (46) 

where  R  is  the  energy  dissipation  rate  required  for  crack 
growth  and  G  is  the  dynamic  energy  release  rate  given  by, 

G  =  if  dU_  dT] 
b\da  da  da) 

.  1  I  dU_  dT) 
bi\dt  dt  dt) 


where  U  is  the  strain  energy,  T  is  the  kinetic  energy,  W  is 
the  work  done  on  the  structure  by  external  loads,  a  is  the 
crack  length,  and  b  is  the  plate  thickness  at  the  crack  tip. 

The  dynamic  energy  release  rate  and  the  dynamic  stress 
intensity  factor  are  connected  through  Freund's  formula. 


[  EG 

U-v" 


AU) 


(48) 


where  G  is  the  dynamic  energy  release  rate,  E  is  Young's 
modulus,  p  is  Poisson's  ratio,  and  A  is  a  function  that  is 
dependent  on  crack  speed.  A  is  given  as 
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where  Cj^  and  C2  are  the  longitudinal  and  shear  wave  speeds. 
[Ref.  2] 

The  wave  speeds  are  given  by  [Ref.  16], 


c. 


lAliH 

P 


(50) 


where  X  and  pL  are  Lcime's  constants,  and 


(51) 
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III.  RESULTS  AMD  DISCUSSION 


A.  PROBLEM  STATEMENT 

The  main  objective  is  how  to  best  approximate  the  behavior 
of  a  continuous  propagating  crack  using  a  discrete  finite 
element  model.  The  model  was  composed  of  a  uniform  mesh  of 
eight  noded  isoparametric  elements.  In  the  first  model  the 
nodes  were  stationary.  To  approximate  the  movement  of  the 
crack,  the  nodal  restraints  were  suddenly  released  at  given 
time  intervals  (time  required  for  the  crack  to  travel  a  nodal 
spacing  in  the  continuous  movement)  .  As  will  be  shown,  this 
resulted  in  spurious  oscillation  in  the  energy  terms.  To 
reduce  the  oscillation  due  to  discrete  jumps  of  the  crack 
movement,  as  simulated  by  the  model,  a  moving  node  element  was 
employed  in  the  second  model. 

Each  of  the  analysis  was  based  upon  the  plane -strain 
conditions  of  a  centered  cracked  plate  of  homogeneous 
isotropic  elastic  material.  The  plate  deformed  under  the 
action  of  an  applied  time- independent  (constant)  traction.  The 
applied  load  was  normal  to  the  crack  plane,  mode  I 
configuration.  The  mass  matrix  was  diagonalized  and  the 
numerical  time  integration  was  accomplished  by  two  different 
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forms  of  the  central  difference  method.  A  comparative  study 
was  conducted  of  the  two  forms. 


B .  VERIFICATION  EXAMPLES 

The  developed  finite  element  analysis  code  was  verified  by 
analyzing  a  problem  with  a  known  solution.  The  problem 
considered  is  governed  by  the  wave  equation.  A  long  uniform 
rod  with  one  end  fixed  and  the  other  end  free  was  -xamined. 
The  free  end  was  subjected  to  an  applied  axial  force  which  was 
released  at  time  zero.  Applying  Newton's  second  law  and  for  a 
constant  area  multiply  by  Young's  modulus  (AE)  the  equation  of 
motion  simplified  to  the  wave  equation 


dC^  dx^ 


(51) 


where  c^  =  E/p  and  u  is  the  axial  displacement  of  an  element 
of  the  rod.  By  the  method  of  separation  of  variables,  the 
analytical  solution  is  expressed  as 


uix,  t)  = 


8FL 


TT 


(-1)" 
(2n+l) 2 


sin  (2n^l)7ix^P3J2n^l)7rct 


(52) 


2L 


2L 


As  shown  in  Figure  4,  the  numerical  result  was  reasonably 
accurate.  To  verify  the  accuracy  of  the  second  model  (moving 
node) ,  it  was  compared  to  the  stationary  model  when  node  five 
was  at  the  center  edge  of  the  element.  For  this  case  the 
results  should  be  and  were  the  same. 
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s 


Figure  4 .  Displacement  of  end  of  beam 
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As  another  comparative  check  of  the  accuracy  of  the 
developed  methods,  the  results  of  the  crack  openina 
displacement  (COD)  was  compared  against  the  analytical 
solution  as  derived  by  Broberg  [Ref.  17]  .  To  obtain  the 
mathematical  description  of  a  two-dimensional  crack,  Broberg 
made  the  foll.,wing  assumption  [Ref.  17]  : 

1.  A  two-dimensional  crack  can  propagate  in  one  plane  only. 

2.  The  material  is  homogeneous  and  isotropic  with  regard  to 
fracturing  characteristics  and  stress- strain  relations. 

3.  The  material  obeys  Hooke's  law  and  is  perfectly  brittle. 

4.  The  surface  energy  is  zero. 

5.  The  crack  propagates  with  a  maximum  velocity  from  the 
start . 

Broberg' s  analytical  solution  applies  to  an  infinite  medium 
with  the  above  characteristics.  These  assumptions  were 
incorporated  into  this  study.  The  comparison  between 
analytical  and  numerical  solutions  are  shown  later. 

C.  NUMERICAL  TIME  INTEGRATION 

An  explicit  procedure  based  on  the  central  difference 
method  was  incorporated  into  this  study.  The  two  forms  of  the 
method  are  as  outlined  in  Chapter  II.  The  latter  is  the 
"summed  form"  [Ref.  13].  From  a  computational  consideration, 
the  "summed  form"  has  an  advantage.  For  a  damped  structure, 
the  "summed  form"  only  requires  the  mass  matrix  be 
diagonalized  to  simplify  the  computation;  whereas,  the  former 
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requires  the  diagonalization  of  the  mass  and  damping  matrix 
for  simplification.  To  date  no  supporting  argument  or  test 
cases  exist  which  justifies  the  use  of  a  diagonal  damping 
matrix . 

During  the  study  it  was  noted  that  the  first  central 
difference  method  was  very  sensitive  to  computational  round 
off  errors.  As  the  time  increment  decreased,  t.he  central 
difference  did  not  converge  while  running  the  model  in  single 
precision.  The  model  was  ran  at  double  precision  (minimizing 
round  off  error)  to  achieve  convergence.  On  the  other  hand, 
the  "summed  form"  converged  during  single  precision  runs. 
Comparison  of  the  convergent  results  revealed  no  difference 
between  the  two  central  difference  techniques.  The  "summed 
form"  simplifies  the  computation,  and  as  a  result,  it  saves 
time  and  money. 

An  important  consideration  when  using  the  central 
difference  method  is  the  time  step.  Because  the  central 
difference  method  is  conditionally  stable,  the  time  step  must 
be  less  than  the  critical  time  step  to  ensure  stability.  The 
critical  time  step  is  defined  as  the  time  required  for  an 
acoustic  wave  to  transverse  the  smallest  element  of  the  system 
[Ref.  11].  A  time  step  less  than  the  critical  time  step 
guarantees  the  stability  of  the  solution.  Stability  meaning 
that  the  displacement  does  not  grow  without  limit.  As 
indicated  in  [Ref.  11],  the  practical  limit  on  At  is 
approximately  25%  less  than  2/0}^,^.  As  seen  in  this  study. 
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decreasing  the  time  step  size  considerably  below  the 
approximate  formula  for  At^^.,  equation  (42j,  resulted  in  no 
change . 

In  the  second  analysis,  the  distance  the  moving  node  can 
travel  is  limited  by  /i .  (3  is  the  fractional  distance  of  the 
element  length.  It  represents  the  minimum  distance  the  moving 
node  should  be  from  the  corner  node.  Figure  5.  Because 
convergence  had  occurred,  decreasing  At  for  a  constant  beta 
resulted  in  no  change  in  the  solutions.  The  minimum  distance 
from  the  corner  node,  /3,  is  dependent  on  At  for  relatively 
large  time  step  size.  As  shown  in  Figure  6,  decreasing  (3  for 
a  constant  At  introduced  higher  frequencies  with  similar 
magnitudes.  If  the  moving  node  propagated  within  the  minimum 
allowable  distance,  the  system  became  unstable.  Exceptionally 
high  frequencies  and  amplitudes  resulted,  Figure  7. 

As  shown  in  the  Table  1,  further  decreasing  the  time  steps 
resulted  in  a  constant  0. 


Table  1.  Fractional  minimum  distance  of 
the  moving  mid -node  from  the  end  node 


NSTEP 

10 

25 

50 

100 

beta  (/3) 

1/5 

1/10 

1/10 

1/10 

*Nstep  is  the  number  of  time  steps  for 
the  crack  to  advance  one  half  the  element 
length. 
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to  element  length 


D.  STUDY  OK  MASS  AND  STIFFNESS  MATRICES 

As  indicated  in  [Ref.  11],  test  cases  have  shown  that  the 
diagonal  matrix,  using  the  procedure  outlined  in  Chapter  II, 
gives  greater  accuracy  than  a  consistent  mass  matrix.  A 
diagonal  mass  matrix  also  requires  less  storage  space  than  a 
consistent  mass  matrix  and  simplifies  the  computation  of  the 
central  difference  scheme. 

For  the  regular  eight  noded  isoparametric  element,  the 
concentrated  mass  at  the  corner  nodes  are  one -eight  that  of 
the  mid-nodes.  For  the  moving  nodal  case,  the  diagonalized 
mass  matrix  from  the  moving  node  element  caused  unstable 
(oscillatory)  conditions.  Because  of  the  close  proximity  of 
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moving  node  five  to  node  one,  the  diagonal  mass  matrix 
procedure  resulted  in  the  majority  of  the  mass  concentrated  at 
nodes  five  and  one.  As  node  five  moved  across  the  edge  of  the 
element  toward  the  center,  the  nodal  point  mass  approached 
that  of  the  regular  element.  Due  to  the  unstable  conditions 
resulting  from  the  diagonal  mass  procedure  for  the  moving  node 
element,  the  diagonal  mass  matrix  was  calculated  with  the 
moving  node  at  the  center  of  the  element  boundary  edge.  The 
computed  diagonal  mass  matrix  remained  constant  throughout  the 
analysis . 

For  the  eight  noded  regular  element  the  element  nodal 
stiffness  remained  constant.  On  the  other  hand,  for  the  eight 
noded  moving  node  element  the  nodal  stiffness  of  nodes  one, 
two,  and  five  varied  as  node  five  moved  along  the  edge  of  the 
element.  Node  five  cared  the  higher  stiffness.  The  remaining 
nodes  (three,  four,  six,  seven,  eight)  were  equivalent  to  the 
regular  element. 

E.  COMPARISONS  OF  ENERGY  TERMS 

In  order  to  make  a  direct  comparison  between  the  use  of  a 
stationary  element  and  a  transition  element,  the  generated 
results  were  compared  with  Broberg's  problem.  The  finite 
element  breakdown  for  Broberg's  problem  was  as  depicted  by 
Kobayashi  et  al.  [Ref.  5]  and  as  shown  in  Figure  8. 

Figure  9  shows  the  crack  opening  displacement  during  crack 
propagation.  The  two  finite  element  results  are  compared  to 
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trhe  analytical  solution  as  solved  by  Broberg.  Except  for  some 
positions,  the  results  were  relatively  similar.  No  comparative 
conclusion  can  be  drawn  from  the  results.  Although  their  close 
comparison  to  the  analytical  solution  indicates  the 
reasonableness  of  the  model. 

Figure  10  shows  the  crack  opening  displacement  at  x  =  0. 
From  this  result,  it  is  seen  that  the  use  of  the  moving  node 
element  approximates  a  propagating  crack  more  closely  than  the 
stationary  element. 

Figure  11  shows  the  work  on  the  plate  during  crack 
propagation.  There  is  a  similar  trend  between  the  stationary 
element  and  the  moving  node  element  with  the  exception  that 
there  is  a  slight  increase  in  the  work  for  the  moving  node 
element . 

Figure  12  shows  the  calculated  strain  energy  during  crack 
growth.  The  stationary  element  produces  spurious  oscillation. 
The  oscillation  is  caused  by  the  instant  release  of  the  crack 
tip  during  the  process  of  discrete  crack  tip  advances  [Ref  5]  . 
As  shown  in  Figure  12a,  there  was  a  rapid  build  up  in  the 
strain  energy  until  the  node  was  released.  Upon  nodal  release 
there  was  a  rapxd  reduction  in  the  strain  energy  followed  by 
a  rapid  increase  until  the  next  nodal  release.  The  oscillation 
amplitude  of  the  strain  energy  increased  during  crack 
propagation.  Employing  a  moving  node  element  made  the  crack 
tip  movement  more  continuous.  As  a  result,  the  oscillation 
amplitude  of  the  strain  energy  was  much  reduced.  Figure  12b. 
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Figure  13  shows  the  calculated  kinetic  energy  during  crack 
growth.  The  delayed  nodal  release,  caused  by  the  stationary 
element,  resulted  in  a  rapid  decrease  in  the  kinetic  energy 
prior  to  the  release  of  the  node,  then  an  rapid  increase 
followed  by  a  rapid  reduction  until  the  next  nodal  release. 
This  trend  was  a  reversal  of  that  seen  for  strain  energy.  The 
oscillation  amplitude  of  the  kinetic  energy  increased  during 
crack  growth  as  it  did  for  the  strain  energy.  Employing  a 
moving  node  element  reduced  the  oscillation  amplitude  because 
the  moving  node  better  represents  a  moving  crack  tip. 

Although  not  necessary  since  the  crack  had  already  passed, 
if  variable  interpolation  was  used  to  move  the  middle  node 
(node  five)  back  to  the  center  position,  there  was  an  increase 
in  the  oscillation  of  the  results.  The  increase  in  oscillation 
is  probably  due  to  the  abrupt  change  in  the  nodal  stiffness 
when  bringing  node  five  back  to  its  center  location.  As  the 
crack  propagated,  there  was  a  gradual  shift  in  the  stiffness 
of  node  five. 

In  an  elastic  medium  there  are  only  two  types  of  waves 
that  propagate  through  a  solid  that  is  unbounded  [Ref.  16]. 
These  are  dilatational  and  distortional  waves.  If  the  solid 
has  a  free  surface,  Rayleigh  surface  waves  may  also  propagate. 
Rapid  reduction  in  stress  occurred  upon  nodal  release.  Due  to 
this  rapid  reduction,  unloading  waves  were  generated  at  the 
crack  tip.  As  the  wave  propagated  through  the  medium,  strain 
energy  decreased  and  kinetic  energy  increased. 
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If  Che  body  is  bounded,  elastic  waves  reflect  and  refract 
off  Che  boundaries  back  toward  the  crack  tip  as  loading  waves. 
The  loading  wave  expands  the  body;  thereby,  amplifying  the 
strain  caused  by  the  tension  load.  This  result  of  an  increase 
in  strain  energy  is  seen  in  Figure  13a. 
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Figure  6.  Effect  of  varying  beta  on  strain  energy  for 
nstepl  =  50:  (a)  beta  =  1/20,  (b)  beta  =  I/IO,  (c)  beta 
1/5. 
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Figure  8.  Finite  element  mesh  for  Broberg's  problem 
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Figure  9.  Crack  Opening  Displacement  for  (a)  stationary 
element  and  (b)  moving  node  element 
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Figure  10.  Crack  opening  displacement  at  x  =  0 
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Figure  11.  Work  vs  Crack  Length  (a)  stationary  element, 
(b)  moving  node  element  (beta=l/5) 
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Figure  12,  Strain  energy  vs  Crack  Length  for  (a) 
stationary  element  (b)  moving  node  element  {beta=l/5) 
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IV.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  study  focused  on  the  finite  element  modeling  of  a 
rapidly  propagating  crack.  Two  methods  were  empxoyed  in  this 
study.  The  first  procedure,  adopting  stationary  elements, 
resulted  in  spurious  oscillation  of  the  energy  terms.  For  this 
procedure  the  propagating  crack  was  represented  by  discrete 
jumps  at  given  time  intervals.  To  more  closely  approximate  the 
propagating  crack,  a  moving  node  procedure  was  employed  in  the 
model.  As  shown  in  the  previous  chapter,  the  moving  node 
procedure  reduced  the  spurious  oscillation  of  the  energy 
terms.  Unlike  the  moving  mesh  procedure,  no  remeshing  is 
required  for  this  procedure. 

From  the  comparison  of  the  two  procedures  mentioned  above, 
the  following  observations  were  noted: 

1.  The  two  forms  of  the  central  difference  methods  yield 
almost  the  same  results.  As  a  result,  the  "summed  form"  was 
preferred  due  to  its  simplicity  and  efficiency  in 
computational  time. 

2.  The  range  which  the  moving  node,  which  represents  the  crack 
tip,  can  propagate  along  the  edge  of  the  element  was  a 
function  of  the  time  step  size  for  relatively  large  time  step 
sizes.  However,  for  reduced  time  step  sizes,  the  range  was  the 
same . 
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3.  For  the  moving  node  procedure,  diagonalizing  the  mass  at 
each  increment  of  movement  resu.'ted  in  sporadic  results.  When 
node  five  is  close  to  the  corner  nodes  the  majority  of  the 
concentrated  element  mass  was  located  at  node  five  and  the 
respective  corner  node  it  was  close  to. 

This  study  did  not  account  for  the  singularity  in  the 
crack  tip.  To  represent  the  stress  field  near  the  crack  tip 
more  accurately,  it  is  recommended  that  singular  element,  as 
derived  in  [Ref.  9],  be  incorporated  into  the  model. 
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